Load data

extractorRData <- function(file, object) {
      #' Function for extracting an object from a .RData file created by R's save() command
      #' Inputs: RData file, object name
      E <- new.env()
      load(file=file, envir=E)
      return(get(object, envir=E, inherits=F))
}

# Load genotypes and phenotypes
load(paste0(data.dir,"ROSMAP_SV_4_association.RData"))

pheno_list_files = list.files(paste0(data.dir,"/",studyDataID,"/association_results"), pattern = "*.assoc.RData", full.names = T)
pheno_list = pheno_list_files[gsub("(.*)/(.*)\\.assoc\\.RData","\\2",pheno_list_files) %in% names(pheno_list)]

res.all = data.frame()
for (pheno_i in pheno_list){
  #pheno_i = pheno_list[1]
  print(pheno_i)
  res = extractorRData(pheno_i, "assoc")
  res_i = res[,c("rs.id","beta","SE", "pval")]
  res_i$bonf = p.adjust(res_i$pval, method = "bonferroni")
  res_i$fdr = p.adjust(res_i$pval, method = "fdr")
  colnames(res_i) = c("sv_id","Estimate","Std. Error", "nom_p", "bonf" , "fdr")
  res_i$pheno = gsub("(.*)\\/(.*)\\.assoc\\.RData","\\2",pheno_i)
  res.all <- bind_rows(res.all,res_i)
}

vcf.meta$sv_info = paste0(vcf.meta$SVTYPE, " chr", vcf.meta$CHROM, ":", vcf.meta$POS, " (len:", vcf.meta$SVLEN, " maf:", vcf.meta$MAF,")")

res_final = vcf.meta[,c("ID","sv_info","closest_gene")] %>% right_join(res.all, by = c("ID"="sv_id")) %>% dplyr::rename("pval" = "nom_p") %>% arrange(fdr)

save(res_final, file = paste0(data.dir,"/",studyDataID,"/association_results.RData"))

Association Results

# Load genotypes and phenotypes
load(paste0(data.dir,"ROSMAP_SV_4_association.RData"))
load(paste0(data.dir,"/",studyDataID,"/association_results.RData"))
phenotypes = phenotypes %>% 
  filter(study==studyID) 

res_final = res_final %>% arrange(pval) %>% 
           filter(ID %in% valid_MAF_SVs) %>%
           filter(ID %in% valid_HWE_SVs) %>%
  filter(pheno %in% names(pheno_list))

res_final = add_LD_info(res_final)

data.table::fwrite(res_final, file = paste0(data.dir,"/",studyDataID,"/association_results.tsv.gz"))
createDT(res_final %>% head(1000))
ci = 0.95
df_m = res_final %>% 
  group_by(pheno) %>%
  mutate(observed = -log10(sort(pval)),
         expected = -log10(ppoints(n())),
         clower   = -log10(qbeta(p = (1 - ci) / 2, shape1 = 1:n(), shape2 = n():1)),
         cupper   = -log10(qbeta(p = (1 + ci) / 2, shape1 = 1:n(), shape2 = n():1)),
         lambda   = median(qchisq(1 - pval, 1)) / qchisq(0.5, 1)) %>%
  reframe(lambda = unique(lambda), lambda_label = sprintf("%s (λ = %.2f)", pheno, lambda)) %>% distinct()

data.table::fwrite(df_m, file = paste0(data.dir,"/",studyDataID,"/lambdas.tsv.gz"))
createDT(df_m)
plot_multitrait_manhattan(res_final, pheno_label = "description")

Session info

sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Stream 8
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.15.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ggrepel_0.9.5     data.table_1.15.4 vcfR_1.15.0       lubridate_1.9.3  
##  [5] forcats_1.0.0     stringr_1.5.1     dplyr_1.1.4       purrr_1.0.2      
##  [9] readr_2.1.5       tidyr_1.3.1       tibble_3.2.1      ggplot2_3.5.1    
## [13] tidyverse_2.0.0  
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.13        ape_5.8            lattice_0.20-45    digest_0.6.36     
##  [5] utf8_1.2.4         R6_2.5.1           evaluate_0.24.0    highr_0.11        
##  [9] pillar_1.9.0       rlang_1.1.4        rstudioapi_0.16.0  vegan_2.6-6.1     
## [13] jquerylib_0.1.4    Matrix_1.6-5       DT_0.33            rmarkdown_2.27    
## [17] labeling_0.4.3     splines_4.1.2      pinfsc50_1.3.0     htmlwidgets_1.6.4 
## [21] munsell_0.5.1      compiler_4.1.2     xfun_0.46          pkgconfig_2.0.3   
## [25] mgcv_1.9-1         htmltools_0.5.8.1  tidyselect_1.2.1   fansi_1.0.6       
## [29] permute_0.9-7      viridisLite_0.4.2  tzdb_0.4.0         withr_3.0.0       
## [33] MASS_7.3-60.0.1    grid_4.1.2         nlme_3.1-153       jsonlite_1.8.8    
## [37] gtable_0.3.5       lifecycle_1.0.4    magrittr_2.0.3     scales_1.3.0      
## [41] cli_3.6.3          stringi_1.8.4      cachem_1.1.0       farver_2.1.2      
## [45] bslib_0.8.0        generics_0.1.3     vctrs_0.6.5        ggeasy_0.1.4      
## [49] RColorBrewer_1.1-3 tools_4.1.2        glue_1.7.0         hms_1.1.3         
## [53] crosstalk_1.2.1    parallel_4.1.2     fastmap_1.2.0      yaml_2.3.10       
## [57] timechange_0.3.0   colorspace_2.1-1   cluster_2.1.2      knitr_1.48        
## [61] sass_0.4.9